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ABSTRACT 


This paper presents a new numerical method for computing incompressible 
turbulent flows. The method is tested by calculating laminar recirculating flows 
fend is applied in conjunction \vith a modified k-e model to compute the flow over 
fe backward-facing step. In the laminar regime, the computational results are 
in good agreement with the experimental data. The turbulent-flow study shows 
that the reattachment length is underpredictfed by,tbe standard k-e model. The 
feddition of a term to the standard model that accounts for the effects of rotation 
On turbulence improves the results in the recirculation region and increases the 
Computed reattachment length. 


Nomenclature 

Aa; viscous operator in the A-directiOn 

Ay viscous Operator in the ^‘•direction 


C. 


pressure coefficient 
constant in the viscosity model (=0.09) 

Cl, C 2 constants in the €-equation(=1.44, 1.92) 

Cs constant for the rotation term in the t-equatiOn 
convection and viscous operator in the x-direction 
convection and viscous operator in the y-direction 
flow quantity 
step height 

block-diagonal matrix, defined in appendix 
=—S{uiUj)/8xj, differenced convective term 
identity matrix 


•Da: 

Dy 

f 

h 

H 

Hi 

I 
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k = turbulence kinetic energy 

Lj length of separation bubble 

i = dissipation length scale 

Nx number of grid points in the a:*directiOn 

Ny number of grid points in the y-direction 

p =h//? + P 

p pressure 

Pj^ production term in A:-eqnation 

=Cie/kPk, production term in e-equation 
Re =tJh/u, Reynolds number 

S source term, defined in appendix 

Sij —^{duifdxj 4- ditjidxi), strain-rate tensor 

Ui mean velocity component in the ^direction 
m'j velocity fiuctuation in the i-directiofi 

U maximum inlet velocity 

Xr reattachment length on Step ■wall 

X4 separation point on no-step wall 

3:5 reattachment point on no-step wall 

d/dxi partial derivative operator 

S/6xi partial differencing operator on a staggered grid 
Ak — — A:”, time-increment of k 

Aq increment vector 

At time-step — f^) 

Ax grid spacing in the x-direction 

Ae — — e”, time-increment of 6 

€ dissipation rate of k 

u laminar kinematic viscosity 

i/g =1 -|“ ut, effective kinematic viscosity 

ut turbulence kinematic viscosity 

p density 

d/j, ( 7 e Pra ndtl n umbers for k and e (= 1 . 0 , 1 . 3 ) 
=z~utu'j 4- %kdij, Reynolds stress 
$ =p -|- O(^), scalar related to pressure 

# shifted cosine transform of ^ 

n rotation term 

Vlij ='^{duifdxj — dUjfdxi), rotation-rate tensor 

Superscript 

* value evaluated at intermediate state 

” value evaluated at time- step Ati 
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Subscript 

i,j indicate direction 

x,y indicate direction a;, y 

1. INTRODUCTION 

Recirculation bubbles are recurring problems in many flows of practical inter- 
est. Tbe flow over a backward-facing step is an example of such flows. In this flow, 
the location of the separation point is fixed and the region where thfe flo\v reattaches 
can be isolated. One can then study the separation-reattachmeht process without 
ahy complexities resulting from motion of the separation point. It is this feature of 
the flow Combined with the simplicity of its geometry that make it a prime candidate 
for both experimental and numerical investigations. 

At the 1980-81 AFOSR-HTTM Stanford Conference on Complex Turbulent 
Flows, 11 groups using 15 methods computed the turbulent flow over a 2:3 sudden 
expansion (Eaton 1982). It was found that all the methods Using the standard 
k-€ model underpredicted the reattachment length as measured by the experiment. 
These results indicate that a modification to the model is necessary. To test different 
models, an accurate numerical method should be used where no artificial viscosity 
is necessary to stabilize the solution. 

One of the objectives of this paper is to present such a method for the incom- 
pressible Navier-Stokes equations. This method is used to test the effect of adding a 
rotation term to the k-e turbulence model on predicting the flow Over the 2:3 sadden 
expansion of Kim, Kline, and Johnston (1980). The particular method (developed 
in §3) uses central differencing on a staggered grid (Harlow and Welch 1965) in 
space and a partially implicit time- advancement algorithm combined with a direct 
Poisson solver to obtain a divergence-free velocity field at each time-step. Results 
for laminar and turbulent computations are presented in §4. 

2. GOVERNING EQUATIONS 


The Reynolds-averaged incompressible Navier-Stokes and continuity equations 


are 


= __p + 


Re dx, 


{Tij -j- 2.^*,) 


( 1 ) 


d 

dxi 


Uj = 0 


( 2 ) 


Here, all the variables are nondimensionalized using the step height h and the 
maximum inlet velocity. All quantities not specifically defined in the text are defined 
in the nomenclature. 
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2.1 Turbulence Mod^l 


To close the above set of equations we need to express the Reynolds stresses, 
Tij, in terms of mean flo’vt' quantities. An eddy-Viscosity model is used: 

‘lU'J'Sij (3) 


■vThere, 




( 4 ) 


Here k is the turbulence kinetic energy and e is the dissipation rate of k. The 
governing equations for k and 6 can be derived from the equations of motion; they 
contain several terms that need to be modeled. The abovfe eddy- viscosity model fOr 
Tij combined with the modeled equations for k and e hrfe known as the k-e model 
(Jones and Launder 1972) . The high-ReynoldS'-nUmber form of the k-e model is 
known to produce poor results in the presence of rotation. Several modifications to 
the model are proposed in the literature to take into account the effects of rotation 
(Rodi 1979; Launder, Priddin, and Sharma 1977; Bardina, Perziger, and Rogallo 
1983). In this work, the term proposed by Bardina, Ferziger, and Rogallo(1983) to 
account for the effects of rotation is Used in conjunction with the k-e equations: 


d 

dt 


^k + -^Ujk 

dt dxj ^ 




Pk — €-i- 


1 d 


Vfjfc dxj j 


Ri dxj V fJfc dxj 




k 


CsOe -f* 


1 d 


/ ilAA 

Re dXj\<T^ dXj ) 


( 5 ) 

( 6 ) 


When 63 = 0 is used, we refer to the model as the standard k-e model. When 
Cs^O we refer to the model as the modified k-e model. 


2.2 Boundary Conditions 

To solve for the system of equations ( 1 - 6 ), boundary conditions are specified for 
all variables except pressure. With second-ordef differencing on a staggered grid, 
the continuity equation at the interior cells, together with the momentum equations 
at the interior grid points and the velocity boundary conditions, provides a closed 
system of algebraic equations for pressure (Moin 1982). 

2 . 2.1 Inlet . All variables except pressure are prescribed at the inlet section. For 
the laminar cases, the streamwise velocity profile {u{y)) is taken to be parabolic. For 
the turbulent cases, the experimental profile (Kim, Kline, and Johnston 1980) for 
u{y) at x/h == 0 is used. The k profile was set to k — Z{u'^ using the 

experimental values of and at xjh = —1.33. The inlet length scale was set 
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(Launder 1982) equal to £ = min{2.5y,0.55}, where y is the distance to the -wall 
and 3 the boundary-layer thickness (e= 0.25h). 

2.2.2 Exit . All Variables except pressure are extrapolated to the exit plane 
using df jdx = 0. 

2.2.3 Wall . At solid-wall boundaries, no-slip is used for the laminar cases. Tur 
turbulent cases, resolution restrictions, as well as the use of the high-Reyndlds- 
number form of the A:*e model, require ut to use near-wall submodels to bridge the 
gap between the wall and the first grid point away from the wall. In this mddel, 
the normal velocity component was set to zero at the wall. To specify the shear 
stress at the wall, we assume that the tangential velocity profile between the first 
grid point away from the wall and the wall is given by the law-of-the-wall (Kays 
1966) as follows: 

u~^ —y~^ y"^ •‘CS 

u-^ =-3.05 + 5/ny+ 5 < <30 

=5.5 -t- 2.5/ny+ 30< y^ <(2000) 

where = u/u*, ahd y~^ = yu*Re. Then given u and y, we can find the value 
of using the Newton iteration method. From the definition of u* we have 

__ 1 du 

p ~ Re dy ~ 


The Neumann boundary condition, dkfdn — 0 was used for the /c-equation. 
The dissipation equation at the cell center adjacent to the wall (point 2) was replaced 
by (Morel et ai. 1981) 


e. = 



(V 


where the dissipation length scale £2 was obtained by interpolation, assuming that 
the dissipation length scale varies linearly between its value on the wall (£1 =* 0) 
and the value at the second cell away from the wall (£3 = 

3. NUMERICAL METHOD 

The equations of motion are discretized on a staggered grid uniform in the x- 
direction but not necessarily uniform in the y-direction. Central differencing is used 
to approximate the spatial derivatives. The velocities are defined on cell surfaces, 
while p, k and e are defined at cell centers. The equations are then advanced in 
time using a fractional step method (Chorin 1968, Temam 1979) combined with the 
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approximate-factoriizatioii technique (Douglas and Gunn 1964; Beam and Warming 
1976; Briley and McDonald 1977). 

8.1 Fractional Step 

To advance the velocity field Ui from time-^ttep il to n-j-l, the following time- 
discretization Of the govefning equations is used: 



At 


^ - if,"-') + 


1 ^ i tif ^ * I 

2Re6xj ^\5xj 



( 8 ) 


u 


n-l-l 




At 


6x. 


■u 


n-+-l 


6 Xi 

= 0 




(&) 

(10) 


Note that the nonlineat terms are advanced by the tecoild-order explicit AdamS- 
Bashforth scheme, "Whereas the diffusion terms are advanced by the second-order 
Crank-Nicholson implicit scheme. Implicit treatment of the viscous terms eliminates 
the numerical viscous stability restriction. The inversion of eqn.(8) would require 
0{N^N^) operations (for — Ny), and OiN^N^) “words of memory; this would 
be costly. This problem is avoided by Using an approximate-factorization method. 
The equations are first Written in A-fotm: 


(/ _ f - <) » - 3r^) + (11) 

The left-hand side of (11) is then approximated as follows! 

(/- fAM-K) - Al(i(3ff? ~ar^)+ (12) 

Note that eqn»(12) is an 0(AJ®) approximation to eqn.(ll). One proceeds by first 
solving a set of tridiagonal equations for w* = (/ -- ^Ay){ul — wf ), and then for 
{u* — uf). The solution of this system of equations requires 0{NxNy) operations, 
a significant reduction in cost. 

3.2 Poisson Equation Solver 

To update the velocity («”"*“ ^) using eqn.(9), has to be determined. This 
is accomplished by applying the numerical divergence operator used in eqn.(lO) to 
6qn.(9): 


Poisson equations such as (13) ai-e efficiently solved using transform methods. We 
first let (Williams 1969) 



for i — 2,3,...,Nx, j = 2,3,...,Ny. This cosine transformation -will enforce a 
boundary condition fof # consistent with the Incompressibility condition alid the 
velocity boundary conditions on daggered grids. Substituting eqh.(14) into eqn.(l3) 
and using the orthogonality property of Cosine, we obtain 



il 

6y^ 




At 6xi 


(15) 


Here k\ = 2[l — cos(j^^^^)]/Aa:^ is the modified wave number, and Ax is the 
uniform grid spacing in the a;-direction. The above tridiagohal system of equations 

can be inverted directly in the y-direction to yield Equation (14) is then used 

to compute The velocity is updated using eqn.(9). Note that p = #-(-0(^), 

in the absence of variable eddy viscosity, i.e =constant, direct substitution yield 

3.3 Update of Viscosity 

In eqn.(8) the effective viscosity (i^e) was lagged so that the k-e equations were 
advanced separately. To integrate the k-e equations in time, the approximate- 
factorization schetne in A-form was used. First the k-e equations are approximated 
using a fully implicit Euler step: 


Ak 

At 


^ (/? "b ) ^n-4-1 


dxj ^ 




Z/j.- 


Reok 6xj 8xj 


n+l^n+1 I 


(16) 


8_^ 


6xj ^ ' ' Recr^dx 




6x, 


(17) 

where Ak = — /c”, and Ae = 6^+^ — e”. The above equations are then 

linearized in time using Taylor series expansion about time-level n to give 


(H AtDx ”i“ AtDy)Aci — AiS 


where 


Aq- 



(18) 

(19) 
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iThe block-diagonal matrix H, and the block vector S ar^ given in the appendix. 
The Dfc and Dy terms are the convection-diffusion operators on a staggered grid in 
the x- and y-direction, respectively. For any reasonable ntimber of grid points the 
inversion of eqn.(18) is prohibitively expensive. This problem is circumvented using 
factorisation to approximate eqn.(16) as follows: 

(H -I- AtD^)n^\E H- AtDy)Aq AtS (20) 

Tquation(20) is an O(Ai^) approximation to eqn.(18). T-^o block tridiagonal matrix 
inversions and a matrix- vector multiplication v^ill yield the solution to eqn^(20). 
Finally, k and £ are updated and the eddy viscosity is computed using eqn.(4). 

4. RESULTS 

The method described in §S was used to simulate tbe flows over bact^ard- 
facing steps in both laminar and turbulent regimes. The geometries of interest are 
those in which the channel walls are parallel so that the apparatus acts as a flow 
diffuser with strong adverse pressure gradient in the Streamwise direction. 

4.1 Laminar flow results. The laminar flow studied experimentally and numeri- 
cally by Armaly et ai. (1983) was chosen for this study. The particular geometry is 
a 1:2 sudden expansion with the inlet channel long enough so that the axial velocity 
profile at the step is parabolic. Computations were carried out using a 130 X 130 
uniform grid for the Reynolds number (He) range 100-800. The present results for 
the reattachment length^ Xr, as a function of Re together with the experimen- 
tal and numerical results of Armaly et ai. are shown in Fig. 1. At He — 600 
the present results start to deviate from the experimental measurements. At this 
Reynolds number, mesh-refinement studies, as well as a check on the effect of the 
location of the exit boundary on the reattachment length, were carried out. From 
these investigations we conclude that the deviation from the experimental results is 
not due to numerical accuracy. A possible explanation for the deviation is that the 
experimental flow becomes three-dimensional at this Reynolds number (as Armaly 
et ai. pointed out). Comparison with the numerical results of Armaly et ai. shows 
that the present results, in the Reynolds number range 600-800, yield a much higher 
feattaohment length than theirs. 

Figure 2 shows the streamlines for Re = 600. Note the appearance of a 
recirculation bubble on the wall opposite the step. Unfortunately, no measurements 
me reported for Re < 1000 for this recirculation bubble. At Re — 1000, the first 
station where measurements are reported, the length of the bubble, L&(= X 4 ), 


^Defined as the point where the velocity changes sign. 
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Figure 1. Reattachment length as a function of Reynolds number in the laminar 
range. 



Figure 2. Streamlines at Re — 600: expansion ratio = 1:2. 


is 10.4 step heights {h). For Re > 1000, the length of the experimentally observed 
bubble decreases with increasing Re. Our computations show that at Re — 600 the 
flow flrst separates on the no-step wall at X 4 = S.5h, and reattaches at 0:5 = 16.3/i, 
resulting in a bubble length = l.Sh. At Re — 800, we compute a bubble growth 
to L 5 = 11.5/i. 
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(a) Case 1, C 3 == 0. 



x/h x/h 

(b) Case 2, C 3 — 0.075. (c) Case 3, C 3 == 0.15. 


Figure 3. Pressure distribution on the no-step side wall; symbols are the experimen- 
tal data (Kim et a/. 1980). 


4.2 Turbulent Flow Results . At high Reynolds numbers, the flow over a backward- 
facing step becomes turbulent, and the mean flow in this regime exhibits one recir- 
culation bubble behind the step. The turbulent flow selected for this study is 
the same as one of the flows selected for thfe 1980-81 AFOSR-HTTM Stanford 
Conference. The geometry is a 2:3 sudden-e^Jpansion (Kim, Kline, and Johnston 
1980), and the Reynolds number at which extensive measurements are provided is 
Re — 44, 580. When this flow was simulated by different groups using the stan- 
dard k-e model, it was found that the reattaehment length is underpredicted by 
the model. Indeed, using the standard k-e model and the boundary conditions 
prescribed in §2, we compute a reattachment length Xr = 5.2/^. The reattachment 
length measured experimentally is about Xr 7/i, with an uncertainty of one step 
height. When we compare our numerical results for the pressure rise on the no-step 
wall with the experimental data (Fig. 3a), we find that a short reattachment length 
results in a shifted (with respect to the experimental data) pressure coefficient, Cp, 
profile. 

These results suggest that a modification to the standard model is necessary. 

It is known that an additional term is needed in the eequation to account for 
the effects of rotation on the length scale (Launder et ai. 1977). It is then not 
surprising that the large recirculation bubble present in the backward-facing step 
flows is mispredicted by the model. If we use the modified k-e model (described 
in §2) with the constant recommended by Bardina et ai. (1983) {C^ = 0.15), 
we find that we overpredict the reattachment length, Xr ~ 9.5/7. By increasing 
the reattachment length, the Cp profile shifts in the desired direction (Fig. 3c). 
However, by overpredicting Xr, the shift is too severe. One can then adjust Cs to 
match closely the experimental reattachment length. We choose to use a constant 
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equal to half the value recommended by Bardina et ai. Fine tuning for a closer 
match was not necessary for the following reasons. Using C 3 = 0.075, the calculated 
reattachment length is Xr = Q.Qh, resulting in good agreement with the data. 
However, if we are interested in the pressure rise on the step-side Wall, we find that 
in all three cases, the calculated results underpredict the experimental curve in the 
range x = 8-16/i (Fig. 4). Fine tuning C 3 will not improve the results in this range. 

In what follows we will refer to cases 1 , 2, and 3 to indicate results obtained 
using C 2 — Oi, 0,075, and 0.15, respectively. With case 1 , the effects of rotation 
as modeled by Bardina et aL are not included. Close examination of Fig. 4 shows 
that accounting for the effects of rotation, case 3, will result in good agreement with 
the experimental data in the range x = 0-7 h. On the other hand, the recovery 
region (a: > 7/t) is poorly predicted. The reasons for the discrepancies can be better 
understood by comparing the calculated mean profiles with the experimental data. 

Figures 5-7 show the mean- velocity profiles at a; = 5.3h, 10.7 h, and 16h for the 
three cases. Clearly, in the recirculation region ( 2 : < 7h, Fig. 5) the maximum 
reverse velocity is best predicted with case 3, and is poorly predicted when the 
rotation term Is not included (case 1 ). Case 2 yields the best overall agreement with 



x/h 


Figure 4. Pressure distribution on the step side wall Symbols are the experimental 
data (Kim et al. 1980), 
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U/Uo VUo U/U„ 


(a) Case 1, C3 = 0. 


(b) Case 2, Cs = 0.075. 


(c) Case 3, C3 = 0.15. /s. 


Figure 5. Mean-velocity profiles at 2 ; = 5.3/i; symbols are the experimental data 
(Kim et a,L 1980). 


the experimental data at a; = 5.3h. The recovery of the mean profile (from one 
showing reverse fiow to a fully developed channel flow profile) is not well predicted 
in all cases (Fig. 6 and 7). In all cases, the calculated recovery is slower than what 



(a) Case 1, C3 = 0. (b) Case 2, Cj = 0.075. (c) Case 3, C3 = 0.15. 

Figure 6. Mean- velocity profiles at a: = lOJh; symbols are the experimental data 
(Kim et al. 1980). 
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(a) Case 1, C3 = 0. 


(b) Case 2, C3 = 0.075. 


(c) Case 3, C3 = 0.15. 


Figure 7. Mean- velocity profiles at x — IQh. Symbols are the experimental data 
(Kim et ai. 1980). 


the experimental data indicate. The good agreelnent observed for case 1 at a: = 
10 Jh (Fig. 6a) is superficial. In case 1, a short reattachment length is calculated, 
so the recovery length at a; = 10.7/i, Lr = 5.5/i(±= x — Xr), over -which the profile 
was adjusting, is much longer than the recovery length over which the experimental 
profile adjusted (L^ = 3.6/i). In other words, if we had used the reattachment point 
as the reference point, the mean profile would not show good agreement with the 
data. 

The profiles of the turbulent kinetic energy and shear stress at three axial 
locations in the recovery region are shown with the experimental data for case 2 in 
Figs. 8 and 9. In general, the location of the peaks and the shape of the profiles 
are well predicted. However, the magnitudes are underpredicted, suggesting that 
the k-e model underpredicts the magnitude of the dissipation length scale in the 
recovery region. 

5. SUMMARY AND CONCLUSIONS 

In this study, a numerical method for computing the incompressible Navier- 
Stokes equations has been developed. The method is time-accurate and the flow 
field satisfies the continuity equation up to machine accuracy at every time- step. 
The method was tested by calculating laminar recirculating flows and was applied 
in conjunction with a k-e model to compute the turbulent flow over a 2:3 sudden 
expansion. 
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Figure 8. Turbulent kintetic energy profiles for case 2; symbols are the experimental 
data (Kim efc ai. 1980). 


In the laminar regime the computational results are in good agreement with the 
data up to Re == 500. For Re > 500, the descrepancy between the computations 
and experiments is due to three-dimensionality in the experimental flow field. The 
turbulent flow study shows that the reattachment length is underpredicted by the 





X = l.lh 


X = 10.3/j 


X = 15.7/1 


Figure 9. Turbulent shear stress profiles for case 2. Symbols are the experimental 
data (Kim et ai 1980). 


14 









standard k-e model. The addition of a term that accounts for the effect of rotation 
to the standard model improves the results in the recircultion region and increases 
the reattachment length. However, the recovery of the mean profiles is still not 
adequately predicted. 

APPENDIX 

The linearization of eqns.(16) and (17) using Taylor series expansions and 
dropping the terms O(At^) and higher yields: 


^ + C„^(2S?+‘Sf/')Ae + At 




dXj ^ 


^ — 


Reak 5xj 
6 


<5 <5 . , 

i/f-t — AA: 


6x. 




Ae 

At 


+ 2C2^Ae - Ca^^A*: + C3n”+'Ae 




k 


' k^ 


1 


ReOt dxi 6x 


VT—Ae-\- - Ca 


+ 


Re(7t dxj 
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6 6 
i/ir— e 




k 


Csn 


(Al) 


n+1. 


6x, 


— 
Sxj ^ 


(A2) 

All quantities without superscript are evaluated at time level n. Equations (Al) and 
(A2) are rearranged in the compact form of eqn.(18) where H is a block diagonal 
matrix with elements 


H, 


( 1 - A(C^2*:(2S?+‘Sf+‘)/e At(C^A^(2S?+'Sr+')/t" + lA 

V-At(C',C„(2Sf+‘Sf+‘)-C2£2/A:2) 1 + At2C2t/A: + 030“+' J 


and S is a block vector with elements 

S* 


^ ^ C^*2(2S?+'S?+')/t - t- 4«,^+‘A + 

lc,C„M2Sf+‘S7+>)- Cat^/A - C3n“+'f- j|;«“+'e+ 


) 


(A4) 
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